RL-TR-95-4 
Final  Technical  Report 
January  1995 


MATSurv  MULTISENSOR  AIR 
TRAFFIC  SURVEILLANCE 


University  of  Connecticut 

Murali  Yeddanapudi,  Yaakov  Bar-Shalom,  Krishna  R.  Pattipati, 
T.  Kirubarajan,  and  Somnath  Deb 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 

19950403  032 


Rome  Laboratory 
Air  Force  Materiel  Command 
Griffiss  Air  Force  Base,  New  York 


This  report  has  been  reviewed  by  the  Rome  Laboratory  Public  Affairs  Office 
(PA)  and  is  releasable  to  the  National  Technical  Information  Service  (NTIS) .  At 
NTIS  it  will  be  releasable  to  the  general  public,  including  foreign  nations. 

RL-TR-95-4  has  been  reviewed  and  is  approved  for  publication. 


APPROVED : 


RICHARD  R.  GASSNER 
Project  Engineer 


FOR  THE  COMMANDER: 


DONALD  W.  HANSON 

Director  of  Surveillance  &  Photonics 


If  your  address  has  changed  or  if  you  wish  to  be  removed  from  the  Rome  Laboratory 
mailing  list,  or  if  the  addressee  is  no  longer  employed  by  your  organization, 
please  notify  RL  (  OCTM  )  Griff iss  AFB  NY  13441.  This  will  assist  us  in  maintaining 
a  current  mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on  a 
specific  document  require  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Pubic  raccrtt-c  burden  for  tHs  cxatecban  d  Wormatfan  is  estimded  togvgraga  1  hod  par  rtsporw,  fr^tring  th.  tine  for 

oah«t^»5™rt^*Ta  th*  cfels  reeded,  ard  corrplethg  trd  reviewing  the  cdectfon  c/  Wormafan  Send  ccmmerts  regadng  the  burden  estmae  a  erry  ahef  asp«*  <* ,his 
Deri.  Mrfiwey,  Sut*  1204,  ArtigaiVA  22202-4302.  end  to  the  Office  Managemert  «rri  Budga  Papawcrit  Redxdon  ProjaX  (070^01  «j.  WasH-gon.  DC  205ta  _ 


AGENCY  USE  ONLY  (Leave  Blank) 


2  REPORT  DATE 

January  1995 


[a  REPORT  TYPE  AND  DATES  COVERED 

Final  May  93  -  Aug  94 


4.  TITLE  AND  SUBTITLE 

MATSurv  MULTISENSOR  AIR  TRAFFIC  SURVEILLANCE 


6.  AUTHOR  (S) 

Murali  Yeddanapudi,  Yaakov  Bar-Shalom, 

Krishna  R.  Pattipati,  T.  Kirubarajan,  and  Somnath  Deb 


5.  FUNDING  NUMBERS 

C  -  F30602-93-C-0090 
PE  -  61102F 
PR  -  2304 
TA  -  E8 
WU  -  P7 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESSES) 

University  of  Connecticut 

Dept  of  Electrical  &  Systems  Engineering 

Storrs  CT  06269-3157 


a  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


N/A 


9.  SPONSORING/MONIT ORING  AGENCY  NAME($)  AND  ADDRESSEES) 

Rome  Laboratory  (0CTM) 

26  Electronic  Pky 
Griffiss  AFB  NY  13441-4514 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


RL-TR-95-4 


1 1 .  SUPPLEMENTARY  NOTES 

Rome  Laboratory  Project  Engineer':  Richard  R.  Gassner/OCTM/ (315)  330-3574 


2a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


12b.  DISTRIBUTION  CODE 


1  a  ABSTRACT  (W«irTrm  200  word.) 

In  this  report,  the  results  of  the  work  on  multisensor  air  traffic  surveillance  and  a 
description  of  MATSurv,  a  software  tool  for  tracking  multiple  targets  using  measure¬ 
ments  from  asynchronous  sensors,  are  presented.  The  work  addresses  the  issues  of 
modeling  the  target  motion,  estimating  target  states,  and  determining  a  consistent 
set  of  measurement-target  associations.  The  specific  objectives  of  this  work  are  to 
formulate  and  solve  the  data  association  problem  for  actual  tracking  scenarios  with 
multiple,  multimodality  sensor  subsystems.  Two  phases  of  the  work  are  described,, 
namely,  the  tuning  phase  whereby  the  measurements  are  associated  with  the  appropriate 
targets  using  the  known  target  IDs  to  tune  and  validate  the  tracking  filter,  and  the 
assignment  phase  whereby  the  association  is  carried  out  using  a  two-dimensional 
assignment  algorithm  developed  at  UConn.  Also,  the  performance  of  the. tracking  and 
assignment  algorithms  using  three  representative  sets  of  measurements  is  illustrated. 
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1  Introduction 


A  surveillance  and  tracking  system  consists  of  a  network  of  geographically  and  functionally  dis¬ 
tributed  sensors  (e.g.,  L-band,  S-band,  C-band,  infrared),  the  objective  of  such  a  system  is  to 
detect  an  unknown  number  of  targets  in  its  field  of  view  and  estimate  the  states  (target  position, 
velocity,  acceleration,  etc.)  using  sensor  measurements  contaminated  by  noise.  This  must  be 
accomplished  in  the  presence  of  spurious  observations  (created  by  background  noise  and  clutter) 
and  occasional  missed  detections  by  sensors. 

The  tracking  process,  as  generally  practiced  today,  consists  of  foui  interrelated  functions. 


1.  Selection  of  the  state  variable  models  used  to  represent  the  target  motion  and  sensor  mea¬ 
surements,  including  models  of  clutter  and  measurements  uncertainties, 

2.  Evaluation  of  an  index  of  desirability  for  each  candidate  measurement-target  association, 

3.  Determination  of  a  consistent  set  of  measurement-target  associations,  and 

4.  Estimation  of  target  states. 


Association  is  the  decision  process  of  linking  measurements  of  a  common  origin  (i.e.,  a  target 
or  false  alarms)  such  that  each  measurement  is  associated  with  only  one  origin.  A  set  of  linked 
measurements  can  then  be  statistically  filtered  to  estimate  the  states  of  targets.  With  the  ever 


increasing  demand  for  higher  performance  in  surveillance  and  tracking  systems,  it  behooves  us  to 


consider  novel  methods  of  associating 
to  estimate  target  states. 


data  from  multiple,  multi-modality  sensor  subsystems,  and 
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As  part  of  the  research  effort  for  Rome  Laboratories,  the  University  of  Connecticut  (UConn)  is 
transitioning  the  data  association  and  estimation  algorithms,  developed  as  part  of  AFOSR  funded 
research,  to  actual  tracking  problems  of  interest  to  the  Air  Force.  The  specific  objectives  of  this 
work  are: 

•  Develop  the  maximum  likelihood  formulation  of  the  data  association  problem  for  actual 
tracking  scenarios  with  multiple,  multi-modality  sensor  subsystems; 

•  Develop  algorithms  for  solving  the  data  association  problem;  and 

•  Demonstrate  the  effectiveness  of  models/ algorithms  on  track  initiation  problem  using  the  RL 
surveillance  testbed. 

We  approached  the  data  association  and  track  initiation  problem  in  two  phases: 

1.  Phase  1:  All  measurements  are  correctly  associated  with  the  appropriate  targets  based  on 
target  ID  (from  the  beacon  returns)  to  validate  the  filter  design  and  fine  tune  the  tracking 
filter  performance.  The  filter  performance  will  provide  a  benchmark  to  evaluate  the  overall 
tracking  and  data  association  algorithm  of  phase  2. 

2.  Phase  2:  Develop  the  maximum  likelihood  formulation  of  the  data  association  problem 
and  solve  the  resulting  problem  using  a  sliding  window  2-dimensional  assignment  algorithm1 
developed  at  UConn.  In  this  phase  the  established  tracks  are  associated  with  the  new 
measurements  from  the  latest  scan  —  this  is  a  2-dimensional  assignment  problem. 

In  this  report,  we  present  the  results  and  a  software  tool  termed  MATSurv  —  Multisensor  Air 
Traffic  Surveillance  System  for  tracking  multiple  targets  using  measurements  from  asynchronous 
sensors. 

JThe  near  optimal  association  obtained  using  the  2D-assignment  algorithm  for  the  measurement  database  pro¬ 
vided,  makes  it  unnecessary  to  use  the  more  general  S-dimensional  algorithm  at  this  juncture.  The  S-dimensional 
algorithm  has  to  be  used  in  more  complex  situations  which  have  crossing,  splitting  and  merging  tracks. 
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2  Description  of  the  Multisensor  Data 

In  this  section  a  description  of  the  raw  scan  data  to  be  used  by  the  sensor  fusion  processor ,  which 

gathers  and  organizes  the  data  from  several  sensors,  is  presented. 

•  The  data  from  the  fusion  processor  consists  of  scans  from  two  L-band  FAA  radars  located 
at  Remsen(“R”),  and  Dansville(“D”),  NY. 

•  The  data  from  two  FAA  radars  consist  of  scans  at  approximately  every  10  seconds.  Each  of 
these  scans  contains  a  number  of  primary  radar  or  skin  returns.  Each  of  these  skin  returns 
consists  of  a  time  stamp,  a  slant  range  and  azimuth  angle  measurements.  For  cooperative 
targets  a  secondary  or  beacon  return  is  also  obtained,  which  provides,  in  addition  to  the 
above,  a  target  identification  number  (squawk)  and  a  target  altitude  measurement. 

•  The  observability  of  the  target  state  requires  a  full  measurement  of  its  position.  Only  beacon 
returns  provide  such  a  measurement  of  the  full  target  position.  All  the  skin  returns  provide 
only  a  partial  measurement  of  the  target  state,  and  hence  are  not  used  in  the  tracking  filter 
at  the  present  stage. 

The  beacon  returns  also  provide  the  target  ID.  We  make  use  of  this  information  as  follows: 

1.  In  the  first  phase  of  the  project,  all  the  measurements  are  associated  with  the  appropriate 
targets  based  on  the  target  ID,  and  the  performance  of  the  tracking  filter  is  evaluated.  This 
phase  helps  to  validate  and  fine  tune  the  filter  performance. 

2.  The  availability  of  the  ID  of  the  measurements  provides  a  means  for  evaluating  the  perfor¬ 
mance  of  the  overall  algorithm  (which  involves  both  association  and  tracking)  in  phase  2. 

Figure  1  shows  the  entire  data  set  of  measurements  with  ID,  available  from  the  two  FAA  radars. 
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3  The  Sensor  Measurement  Statistics 


The  measurement  noises  of  the  sensors  are  assumed  to  be  zero-mean  white  with  variances  as 
follows. 

For  the  range,  the  resolution  cells  are 

ct  ,  . 

Ar  =  y,  (1) 

where  c  =  3  ■  108  m/s  and  r  is  the  pulse  width  (6  [is  for  Remsen  and  1.8  [is  for  Dansville).  It  is  as¬ 
sumed  that  the  range  measurement  is  uniformly  distributed  in  each  resolution  cell,  hence  the  stan¬ 
dard  deviation  of  the  range  measurement  noise  is  ar  —  ^|,  which  yields  0.5196  km  =  0.2806  nmi 
for  Remsen  and  0.1559  km  =  0.0842  nm.i  for  Dansville. 

The  azimuth  measurement  noise  standard  deviations  were  taken  as  ag  =  2.618  mrad  =  0.15° 
based  on  FAA  data.  No  other  specific  information  was  available. 

The  altitude  measurement  noise  was  taken,  based  on  FAA  data,  as 


4  Design  of  the  Tracking  Filter 


A  description  of  the  design  and  implementation  details  of  the  tracking  filter  is  given  below.  The 
overall  block  diagram  of  this  filter  is  shown  in  Figure  2. 


Figure  2:  Block  diagram  of  the  tracking  filter  with  conversions  between  three  coordinate  systems 

(sensor,  target-local  and  the  global). 


The  following  steps  in  the  tracking  filter  design  are  described: 

•  The  scenario  of  interest  cannot  be  modeled  as  linear  both  due  to  the  nonlinearity  of  mea¬ 
surement  equations  as  well  as  to  the  large  geographical  area  involved.  The  latter  makes  it 
necessary  to  use  a  spherical  earth  model  in  place  of  a  flat  earth  approximation. 

•  In  order  to  be  able  to  apply  the  standard  Kalman  filter,  we  need  a  three  dimensional  Carte¬ 
sian  coordinate  frame  of  reference.  In  addition,  the  instantaneous  target  motion  in  this 
frame  of  reference  is  approximated  by  a  linear  model. 
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We  start  the  description  of  the  tracking  filter  after  the  kth  sampling  interval.  The  esti¬ 
mated  target  state  in  geographical  coordinates  xg(k\k)  is  a  8  x  1  vector  with  the  following 
components: 

x^,(k\k)  -  latitude,  [— 7t/2,  7t/2] 

x\(k\k)  -  longitude,  [ — tt ,  7 r] 

xa(k\k)  -  altitude  in  km  (above  MSL) 

x^(k\k)  -  velocity  due  north  in  km/s 
xs(k|&)  =  (2) 

x\(k \k)  -  velocity  due  east  in  km/s 
xa(k\k )  -  vertical  velocity  up  in  km/s 
x^{k\k)  -  acceleration  due  north  in  km/s2 
x\(k\k)  -  acceleration  due  east  in  km/s2 

As  can  be  seen  from  the  above  definition,  the  state  estimate  xg(k\k )  is  in  mixed  dimen¬ 
sions  (i.e.,  in  angles  and  distances).  Nevertheless,  we  define  the  covariance  matrix  Pg(k\k) 
completely  in  terms  of  distance  and  distance  rates  only.  The  compatibility  between  the 
state  in  mixed  dimensions  and  the  covariance  in  uniform  dimensions  can  be  understood  by 
interpreting  latitude  as  distance  due  north  of  equator  and  longitude  as  distance  due  east  of 
Greenwich  meridian.  We  can  thereby  define  Pg(k\k)  as  the  covariance  matrix  associated  with 
the  distance  errors  due  north,  east  and  up  and  errors  in  their  corresponding  time  derivatives.. 

A  “local  Cartesian  coordinate  frame”  is  defined  based  on  the  target  state  estimate  xg(k\k). 
The  origin  of  this  reference  frame  is  [x^,(k\k)  x\(k\k)  0]  and  the  axes  are  oriented  along 
the  north  n^,  east  e^,  and  altitude  directions.  The  symbol  Lk  is  used  to  denote  this  local 
reference  frame  uniquely  defined  by  the  origin  and  the  three  coordinate  axes. 


•  The  raw  measurement  z(k  +  1)  vector  (consisting  of  slant  range,  azimuth  and  altitude)  is 
first  transformed  into  the  geographical  coordinates  z g(k  +  1).  The  three  components  of 
zg(k  +  1)  are  the  latitude,  longitude  and  altitude.  The  raw  measurement  covariance  R(k  +  1) 
is  transformed  to  Rg(k  +  1)  via  the  first  order  approximation:  Rg  =  [VzZ^  R  [vzz^j  .  The 
Jacobian  Yzz  evaluated  at  the  raw  measurement  z(k  +  1). 

•  The  measurement  zg(kpl)  and  covariance  Rg(k+l)  are  transformed  into  the  local  Cartesian 
axes  Lk.  This  transformation  yields  zLk(k  +  1)  and  RLk(k  +  1).  Note  that  the  measurement 
zg(k  +  l)  is  in  mixed  dimensions  (two  angles  and  one  distance)  whereas  zik(k  +  1)  consists  of 
all  distance  components.  A  similar  mapping  occurs  in  the  transformation  of  the  measurement 
covariance  matrix  from  Rg(k  +  1)  to  RLk{k  +  1). 

•  The  geographical  state  estimate  xg(k\k)  and  covariance  Pg(k\k )  are  mapped  into  the  local 
coordinate  axes  Lk.  This  yields  XLk{k\k)  and  Pik{k \k)  as: 

r  ~\ 1 

x Lk{k\k )  =00  xa(k\k)  x^,(k\k)  x\(k\k)  xa(k\k)  x^(k\k)  xx(k\k)  (3) 

PlMk)  =  Pg(k\k)  (4) 

The  predicted  state  x^k{k  +  1| k)  and  covariance  PLk(k  +  l\k)  are  determined  by  assuming 
that  the  instantaneous  target  motion  is  linear  along  the  local  coordinate  axes  Lk- 

xLk(k  +  l\k)  =  $jkXLjt(fc|fc)  (5) 

PLk(k  +  l\k)  =  $kPLk(k  +  l\k)$'k  +  rkQkT'k  (6) 


S 


•  The  choice  of  an(k),  ae(k)  and  au(Jc)  is  a  critical  design  issue.  Since  an(k)  and  ae(k) 
are  process  noise  levels  along  north  and  east  axes,  they  are  set  to  the  same  value,  i.e., 
an(k)  =  ue(k)  =  <j k{k)  where  h  denotes  horizontal.  The  vertical  motion  of  the  target  is 
more  predictable  than  along  the  horizontal  axes.  Hence,  the  process  noise  level  along  the 
altitude  axes  aa(k)  is  different  from  Oh{k ).  Both  the  horizontal,  ah(k ),  and  the  altitude, 
cra(k ),  process  noise  levels  are  chosen  according  to  the  graph  shown  in  Figure  3. 


2The  time  interval  can  in  general  be  negative,  particularly  when  the  detection  at  tk+i  and  the  detection  at 
tk  originate  from  different  radars.  The  implementation  of  the  filter  for  negative  time  updates  is  explained  in  the 
appendix. 
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Figure  3:  Process  Noise  Level  as  a  function  of  Sampling  Interval 


•  For  c Th(k )  we  chose 


<7hm„  =  7.5  m/s2 
Sh.  =  0.0  s 

"'min 

and  for  <7a(k)  the  following  values  were  chosen 


c h  =  0.0  m/s 2 

umm  / 

4max  =  40.0  5 


^&min 


0.5  m/.s2 

1.0  5 


=  0.0  m/s 2 
40.0  s 
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•  The  linear  Kalman  filter  state  and  covariance  update  equations  are  implemented  in  the  local 
Cartesian  coordinates.  These  update  equations  yield  X-ik{k  +  l|fc  +  1)  and  Pik(k  +  l\k  +  1). 
In  addition  we  get  the  innovations  VLk{k  +  1)  and  the  innovations  covariance  SLk(k  +  1) 
which  are  used  to  gauge  the  performance  of  the  filter. 


vLk{k+l) 

=  zLk(k  +  1)  -  H  xLk(k  +  l\k) 

(9) 

S Lk(k  +  1) 

=  H PLk(k  +  l\k)  H' +  RLk(k  +  l) 

(10) 

WLk(k  +  1) 

=  PLk(k  +  l\k)H'SLk(k  +  l)-1 

(11) 

x Lk(k  +  1|  k  +  1) 

=  *Lk(k  +  l\k)  +  WLk(k  +  1)  vLk(k  +  l) 

(12) 

PLk{k  +  l\k  +  1) 

=  pLk(k  + 1|  k)-  wLk(k  + 1)  sLk(k  + 1)  wLk(k  +  iy 

(13) 

where  H  =  [I  0  0]3x8  and  the  8  X  3  matrix  WLk(k  +  1)  is  the  filter  gain. 

•  The  updated  state  x.Lk(k  +  1| k  +  1)  and  covariance  PLk(k  +  1| k  +  1)  are  along  the  Lk 

reference  frame.  To  complete  the  cycle  of  transformations,  we  transform  the  updated  state 
xLfi(k+l\k+l)  and  covariance  PLk(k+l\k+l)  into  geographical  coordinates  xff(/M-l|&+l)  and 
Pg(k  +  l\k  +  1),  respectively.  The  two  equations  below  illustrate  the  sequence  of  operations 
that  are  involved  in  this  transformation. 

xLfc(^’  +  1| k  +  1)  =>  x5(&  +  1| k  +  1)  = =>-  Lk- i-i 

PLk(k  +  l\k  +  l)  Lk. ^+1  xg(k  +  l\k  +  l),  Pg(k  +  l\k  +  l) 

•  The  filter  is  initialized  using  the  first  two  measurements,  i.e.,  xa(2|2)  and  Pfl(2|2)  are  de¬ 
termined  using  two  point  diferencing  from  zg(l),  Rg(  1)  and  zs(2),  Rg( 2).  Since  only  two 
measurements  are  used  the  initial  acceleration  components  x^(2\2)  and  xa(2|2)  are  both  set 
to  zero. 
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5  Data  Association 


In  this  section  we  present  the  formulation  and  results  of  the  data  association  problem.  A  sliding 
window  2-dimensional  assignment  algorithm  has  been  used. 

The  formulation  of  the  2-D  assignment  problem  is  presented  below.  Consider  the  scenario  at 
scan  k  >  1  .  There  are  n(k)  validated  tracks  from  the  previous  assignments,  and  corresponding  to 
track  i,  i  =  1,  •  •  •  ,  rc(h),  we  have  the  latest  state  estimate  x^(&|&)  and  covariance  P^(k\k).  Let 
scan  k  -f  1  contain  m{k  +  1)  measurements. 

The  hypotheses  upon  which  the  2-D  assignment  algorithm  is  based  on  are  the  following: 

1.  Measurement  j ,  1  <  j  <  m(k  +  1)  originated  from  a  target  corresponding  to  one  of  the  n(k) 
validated  tracks,  say  track  i,  i  =  1,  -  •  •  ,n(k). 

2.  Measurement  j,  1  <  j  <  m(k  +  1)  is  a  false  alarm.  Track  index  i  —  0  is  used  to  designate  a 
dummy  target  (i.e. ,  a  source  of  false  alarms).  Such  a  measurement  is  also  kept  as  a  candidate 
for  a  new  track,  i.e.,  if  within  Tq  =  30s  it  has  another  measurement  in  its  neighborhood  it 
initiates  a  new  track. 

Each  track  (excluding  track  0)  is  assigned  at  most  one  measurement,  and  each  measurement  is 
assigned  to  at  most  one  track.  On  the  other  hand,  the  number  of  measurements  that  may  be 
assigned  to  track  0  is  not  limited. 

An  existing  track  is  dropped  if  within  Tjj  =  60s  no  new  measurement  is  associated  with  it. 
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5.1  Feasible  Assignments 

We  shall  define  an  assignment  in,  as  the  mapping  between  the  measurement  indices  j  and  the 
track  indices  i.  This  assignment  can  be  represented  using  the  set  of  binary  (0, 1)  variables  {pi,j}- 
A  feasible  assignment  must  satisfy  the  following  requirements 

w  =  {pi,j  €  {0,1},  *  =  0,...,n(k),  j  =  1,. m(k  +  1)}  (14) 

n(k) 

Y^Pi  j  =  1  j  =  l,...,m(fc  +  1)  (15) 

t=0 

m(/c-fl) 

=  1  i  —  l,  ■  •  ■ ,  n(k)  (16) 

i=x 

5.2  The  2D  Assignment  Problem  and  its  Complexity 

Let  be  the  set  of  all  feasible  assignments,  then  the  total  number  of  feasible  assignments,  i.e. , 
the  cardinality  of  the  set  fi,  can  be  shown  to  be 

M  =  t  (f)  cp,  (17) 

where  p  =  min {n{k),m{k  +  1)}  and  q  =  ma x{n(k),m(k  +  1)}. 

Let  Cij  be  the  cost  incurred  in  assigning  measurement  j  to  track  i.  The  2-dimensional  assign¬ 
ment  problem  can  now  be  stated  simply  as  finding  an  optimal  assignment  u>*  that  minimizes  the 
overall  cost,  i.e., 

n(k) 

w*  =  argmin£  pijaj  (18) 

1=0  j= 1 

The  number  of  feasible  assignments,  |fl|,  is  very  large  even  for  moderate  values  of  p  and  q,  hence  an 
efficient  assignment  algorithm  is  required  for  solving  this  problem.  The  modified  auction  algorithm 
is  ideally  suited  for  solving  this  2-dimensional  assignment  problem.  A  detailed  description  of  this 
algorithm  can  be  found  in  [1], 

In  the  present  context  of  assigning  measurements  to  tracks,  the  cost  criterion  ct  J-  is  the  negative 
logarithm  of  the  ratio  of  the  likelihood  of  measurement  j  originating  from  track  0)  to  the 
likelihood  of  measurement  j  originating  from  track  0  (i.e.,  the  likelihood  that  measurement  j  is  a 
false  alarm). 
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The  evaluation  of  c2j  is  presented  in  the  following.  The  likelihood  that  measurement  j,  1  < 
j  <  m(k  +  1),  originated  from  track  i,  1  <  i  <  n(k),  is  given  by 

A i,j(k  +  1)  =  Pd  ■  |27t S{j(k  +  1)|_5  .  exp  | -~rn,j(k  +  1)}  (19) 

where  P 'p  is  the  probability  of  detection,  rju](k  +  1)  is  the  normalized  innovation  squared  given  by 

Vi,j{k  +  1)  =  ui,j{k  +  1)  [Si,j(k  +  1)]  ,yi,j(k  +  1 )  (20) 


Vij(k  +  1)  is  the  innovation  and  Sij(k  +  1)  the  innovation  covariance,  associated  with  the  track  i 
filter  updated  with  measurement  j. 

We  shall  assume  that  false  alarms  are  uniformly  probable  in  the  whole  of  the  surveillance 
region  of  volume  ']/(&  +  1).  Since  false  alarms  are  assigned  to  track  0,  the  likelihood  that  the 
measurement  j  is  a  false  alarm  is, 


Ao,#  +  1)  — 


(21) 


$(k  +  1) 

The  cost  C{j  of  assigning  measurement  j  to  track  i  is  the  negative  log-likelihood  ratio,  given  by 


A  j,j(k  +  1) 
;Ao  ,j(k  +  1)J 


=  7^Ak  +  I)  +lo§ 


'\2vShj{k  +  1)1§ 

pDy(k  + 1) 


(22) 


The  following  values  have  been  used  for  the  probability  of  detection  Pp  and  the  surveillance  region 
volume  of  the  region  in  which  a  candidate  false  alarm  is  uniformly  distributed: 


PD  =  0.999 

(23) 

T(fc  +  1)  =  1000  km3 

(24) 
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6  Examples 

The  performance  of  the  tracking  filter  and  the  association  algorithm  is  illustrated  for  two  typical 
target  trajectories. 

In  example  1,  the  association  algorithm  perfectly  matches  the  measurement  ID  data,  and  hence 
the  output  of  the  tracking  filter  with  ID  and  with  association  is  quite  identical.  These  results  are 
shown  in  figures  4-11.  Note  that  the  tracking  filter,  handles  the  target  maneuver  (which  is  a 
typical  landing  pattern)  successfully. 

In  example  2  the  trajectory  formed  using  the  association  algorithm  is  definitely  superior  to  the 
one  formed  using  IDs.  The  results  for  this  example  are  shown  in  figures  12-23. 

This  case  is  unusual  because  it  appears  that  the  ID  was  changed  toward  the  end  of  the  flight 
(before  landing).  The  ID-based  association  came  up  with  19  measurements  while  the  assignment 
algorithm  obtained  the  same  19  plus  another  6  by  which  time  there  was  a  new  ID  but  the  mea¬ 
surements  gave  a  good  filter.  This  shows  the  superiority  of  our  assignment  algorithm  over  the 
“ground  truth”  which  has  its  quirks.  Once  again  note  that  the  tracking  filter  is  quite  robust  in 
handling  the  rather  sharp  target  maneuver. 

Example  3  illustrates  two  very  positive  aspects  of  the  association  algorithm.  Firstly  the  asso¬ 
ciation  algorithm  drops  a  measurement  that  is  very  noisy  (i.e.,  outlier  rejection),  thus  yielding  a 
better  track  than  that  is  obtained  using  the  ID.  Secondly,  its  brings  in  a  measurement  that  has 
been  (incorrectly)  assigned  an  ID  corresponding  to  a  different  target.  Comparing  the  tracks  shown 
in  figures  27  and  28,  it  can  be  seen  that  dropping  the  outlying  measurement  definitely  improves 
the  track.  The  normalized  innovations  squared  shown  in  figures  31  and  32  clearly  indicates  that 
the  measurement  that  is  dropped  is  an  outlier  while  the  new  measurement  (with  a  different  ID) 
definitely  belongs  to  this  target.  This  is  a  known  phenomenon  of  transponder  interference  when 
two  aircraft  are  close.  In  this  case  the  “ground  truth  might  not  be  true  (ID  codes  have  been 
switched  due  to  the  interference). 
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MAI  Surv :  Multisensor  Air  Irattic  Surveillance  Software 


Scans  Iracks  Iargets  Sector  Sensor  filter  Beset 


Figure  4:  Example  1.  Target  trajectory,  using  ID 


MATSurv:  Multisensor  Air  Traffic  Surveillance  Software 


£cans  Iracks  Jargets  Rector  Sensor  filter  Beset 

Target  Index:  14  (67) 

Target  IFF  Code:  21 11 
Number  of  Detections  :  31 
:  Average  sampling  tfme  Interval :  8.8645  seconds 

Minimum  time  interval  between  successive  detections  :  0.6992  seconds 

Maximum  time  Interval  between  successive  detections  :  30.3008  seconds 
Time  Interval  between  the  first  artd  thelast  detections :  273.8008  seconds 
1  detections  received  out  of  sequence. 

Longest  backtrack :  1.0000  seconds. 

I :  Updated  Target  Position 
h  :  Predicted  Target  Position 
* :  Measurement  from  Remsen  Radar 
* :  Measurement  from  Dansville  Radar 


Figure  6:  Example  1.  Target  trajectory  information,  using  ID 


MATSurv:  Multisensor  Air  Traffic  Surveillance  Software 


Scans  Tracks  Targets  Sector  Sensor  filter  Beset 


Track  Index  :  27  (231)  ; 

Number  of  Detections :  31 

Average  sampling  time  interval :  8.8000  seconds 

Minimum  time  interval  between  successive  detections  :  0.0000  seconds 

:  Maximum  time  interval  between  successive  detections  :  30.3008  seconds 
Time  interval  between  the  first  and  the  fast  detections :  272.8008  seconds 
1  detections  received  out  of  sequence. 

Longest  backtrack :  -0.0000  seconds. 

The  31  detections  originate  as  follows: 

31  detections  from  target  with  Iff:  211 1 


Figure  7:  Example  1.  Target  trajectory  information,  using  association 


Figure  10:  Example  1.  Output  of  tracking  filter:  Normalized  innovation  squared  (in  3  dimensions) 


Figure  11:  Example  1.  Output  of  tracking  filter:  Negative  log-likelihood  ratio 
(should  be  below  0  to  indicate  “target”  more  likely  than  “false  track”) 


Scans  Tracks  Iargetg  Sector  Censor  mixer  _ . 

Target  Index :  39  (67) :  i  : 

Target  IFF  Code :  08 
Number  of  Detections :  19 

Average  safnpHng  time  Interval :  10.6423  seconds 
Minimum  time  Interval  between  successive  detections  :  9.0000  seconds 

Maximum  time  Interval  between  successive  detections :  20.3000  seconds 
Time  Interval  between  the  firs^  arid  the  last  detections :  202.2031  seconds 
0  detections  received  out  of  sequence. 

Longest  backtrack:  -0.0000  seconds. 

i :  Updated  Target  Position 
* :  Predicted  Target  Position 
* ;  Measuremeni  fromiRemsen;  Radar 
* :  Measurement  from'Dansville  Radar 


Figure  15:  Example  2.  Target  trajectory  information,  using  association 
Note:  25  detections  with  old  &  and  new  ID 


Figure  16:  Example  2.  Estimated  trajectory  of  target,  using  ID 

(19  detections) 


l  lH  MAT  Surv :  Mullfsensor  Air  Traffic  Surveillance  Software  111111 


Figure  17:  Example  2.  Estimated  trajectory  of  target,  using  association 

(25  detections) 
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Figure  19:  Example  2.  Estimated  altitude  of  target,  using  association 

(25  detections) 


Figure  20:  Example  2.  Output  of  tracking  filter:  Normalized  innovation  Squared,  using  ID 

(19  detections) 


;  MAT  Surv :  Multisensor  Air  Traffic  Surveillance  Software  OfSI j 


Figure  21:  Example  2.  Output  of  tracking  filter:  Normalized  Innovation  Squared,  using  association 

(25  detections) 
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MAISurv:  Mulhsensor  Air  Irallic  Surveillance  Sot  I  ware 


Scans  Iracks  Iargets  Sector  Sensor  Fitter  Beset 


100  120  140  160  ISO  200  220 


et  39.  Iff :  88.  with  1 9  detections.  Negative  LoqllkellhOoti  Ratio  vs.  Time. 


Figure  22:  Example  2.  Output  of  tracking  filter:  Negative  log-likelihood  ratio,  using  ID 

(19  detections) 


Figure  23:  Example  2.  Output  of  tracking  filter:  Negative  iog-likelihood  ratio,  using  association 


(25  detections) 
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MATSurv:  Multiscnsor  Air  Traffic  Surveillance  Software 


Scans  Iracks  Targets  Sector  Sensor  Alter  Beset 


'Target  Index :  7  (6?)  • 

Target  IFF  Code  :3019  | 

Number  of  Detections :  40  1 

Average  sampling  time  Interval :  7.1400  seconds 

Minimum  time  Interval  between  successive  detections  :  1.0000  seconds 

Mcbdmum  time  interval  between  successive  detections  :  1 1.9023  seconds 
Time  Interval  between  the  first  and  thelast  detections :  281.1992  seconds 
3  detections  received  out  of  sequence. 

Longest  backtrack :  2.1016  seconds. 

i :  Updated  Target  Position 
* :  ’Predicted  Target  Position 
*  :  Measurement  irom  Remsen  Radar 
* :  Measurement  from  Dansville  Radar 


Figure  26:  Example  3.  Target  trajectory  information,  using  ID 
40  detections,  all  detections  have  the  same  ID  3019 


pi 


MATSurv:  Multisensor  Air  Traffic  Surveillance  Software 


Scans  Tracks  Targets  Sector  Sensor  Filter  Reset 


Track  Index: 9 (231)  1 

Number  of  Detections :  40  | 

I 

Average  sampling  time  Interval :  6.9199  seconds 

Minimum  time  interval  between  successive  detections  :  0.0000  seconds 

|  Mcjidmum  time  Interval  between  successive  detections :  1 1 .9023  seconds 
1  Time  interval  between  the  first  and  the  last  detections :  276.7969  seconds 
3  detections  received  out  of  sequence. 

Longest  backtrack :  -0.0000  seconds. 

The  40  detections  originate  as  follows: 

39  detections  from  target  with  iff:  3019 
1  detections  from  target  with  iff:  1521 


Figure  27:  Example 
40  detections:  one  detection  with 


3.  Target  trajectory  information,  using  association 
ID  3019  is  discarded,  and  one  detection  with  ID  1521  is  added 
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Scans 

4Z*15‘N 

Jracks  Iargeis  Sector  Sensor  filter  Beset 

111 

76*57*W  76-5VW  76’45'W  76*39’W  76’33'W  76*27^  76*21*W  76‘15‘W  76*09'W  76*03'W 

Target  7.  iff :  3019.  with  40  detections.  Filtered  Latitude  vs.  Longitude,  with  error  ellipses. 


Figure  28:  Example  3.  Estimated  trajectory  of  target,  using  ID 
Observe  the  effect  of  the  outlying  detection  in  the  top  right  corner  of  the  track 


Notice  that  the  outlying  detection  (with  ID  3019)  has  been  discarded 
A  new  detection  (with  ID  1521)  has  been  included  in  the  top  right  corner 


Figure  30:  Example  3.  Estimated  altitude  of  target,  using  ID 


Figure  31:  Example  3.  Estimated  altitude  of  target,  using  association 


Note  that  the  additional  detection  (with  ID  1521)  “conforms”  with  the  rest  of  the  target  track 
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Figure  32:  Example  3.  Output  of  tracking  filter:  Normalized  Innovation  Squared,  using  ID 
The  sharp  spike  in  the  NIS  is  due  to  the  presence  of  the  outlying  detection 


7  Discussion 


In  this  section  we  shall  discuss  the  major  results  of  phase  1  and  phase  2  of  this  project.  An 
interactive  software,  MATSurv  :  Multisensor  Air  Traffic  Surveillance,  that  runs  in  MS-Windows 
(version  3.0  or  higher)  has  been  developed  to  analyze  the  performance  of  the  association  algorithm 
and  display  the  results  in  a  graphical  format. 

The  results  of  phase  1  relate  to  the  performance  of  the  tracking  filter  applied  to  measurements 
classified  using  the  ID  (designated  as  if  f  code  in  the  database).  These  results  helped  in  arriving  at 
the  best  choice  of  the  design  parameters  for  optimal  filter  performance.  The  measurement  database 
(in  file  data.bin)  from  the  two  FAA  radars  contained  detections  of  targets  that  were  in  a  variety 
of  trajectories.  While  many  of  these  target  trajectories  could  be  described  by  a  straightforward 
2nd  order  motion  model  with  a  low  process  noise,  there  were  some  maneuvering  targets  that  would 
require  at  least  a  2,ld  order  (or  3nd  order)  motion  model  with  considerable  process  noise.  Since 
the  same  tracking  filter  is  required  to  handle  these  two  extreme  cases,  the  choice  of  the  design 
paramenters  has  to  be  necessarily  conservative,  i.e.,  tuned  to  handle  the  worst  case.  This  requires 
trading  off  some  of  the  achievable  estimation  accuracy  for  an  enhanced  ability  to  track  targets 
during  maneuver. 

The  above  observation  clearly  indicates  that  replacing  the  Kalman  filter  (the  central  block  in 
the  tracking  filter  shown  in  figure  4)  with  an  Interacting  Multiple  Model  filter[2]  would  make  the 
above  design  trade-off  unnecessary  and  will  enhance  the  performance. 

In  phase  2,  the  measurements  were  stripped  of  their  IDs  and  processed  using  our  assignment 
algorithm.  The  results  obtained  indicate  that  the  association  algorithm  provides  a  superior  clas¬ 
sification  of  the  measurements  into  tracks  (i.e.,  trajectories  of  the  hypothesized  targets)  than 
compared  to  the  target  trajectories  obtained  using  the  measurements  IDs.  The  multiplicity  of 
targets  assigned  to  the  same  ID  prevents  the  exclusive  reliance  on  the  target  ID,  and  its  use 
in  evaluating  the  performance  of  the  association  algorithm  is  clearly  inappropriate.  A  particu¬ 
lar  track  formed  by  the  association  algorithm,  has,  in  general,  a  few  measurements  less  than  the 
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corresponding  target  trajectory  obtained  using  the  IDs.  This  is  primarily  due  to  the  fact  that 
the  association  algorithm  rejects  measurements  (i.e.,  outliers )  that  deviate  considerably  from  the 
established  track.  Discarding  these  measurements  yields  a  better  estimate  of  the  trajectory  than 
the  one  obtained  by  including  these  outlying  measurements. 
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8  Future  Work 


The  following  are  some  of  the  major  issues  that  are  worthy  of  attention  in  the  future. 

•  Incorporation  of  partial  2D  measurements  (skin  returns)  into  the  overall  association-estimation 
framework.  These  2D  measurements  can  be  used  in  the  enhancing  the  purity  of  an  estab¬ 
lished  track  formed  using  the  full  3D  measurements  (beacon  returns). 

•  The  use  of  an  Interacting  Multiple  Model  estimator  to  handle  both  maneuvering  and  non¬ 
maneuvering  targets  with  maximum  accuracy. 

•  Real-time  operation. 

•  The  2-dimensional  assignment  algorithm  used  in  the  present  context  can  be  extended  to  the 
more  general  S-dimensional  case.  The  near  optimal  results  obtained  using  2D  assignment 
for  the  currently  available  measurement  database,  suggest  that  in  order  to  fully  utilize  the 
advanced  features  of  the  S-D  assignment,  a  measurement  database  that  contains  a  more 
complex  scenario  is  needed. 

•  Of  more  theoretical  interest  is  the  issue  of  modifying  the  currently  used  ML  formulation  of  the 
data  association  problem  to  a  more  general  MAP  or  MMSE  kind  of  a  formulation.  This  would 
require  further  modifications  to  the  auction  algorithm  wherein  successive  2D  assignments  are 
aggregated  into  MAP  trajectory  estimates  based  on  the  measurement  data. 
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A  Negative-Time  Update 

In  this  section  we  outline  the  negative  time  update  procedure  used  in  our  implementation  of  the 

tracking  filter.  In  this  case  the  time  interval  Sf.  =  tk+i  —  tk  is  negative.  In  the  standard  case  (i.e., 

8k  >  0)  the  past  is  completely  described  by  the  current  estimate  (i.e.,  it  is  the  sufficient  statistic): 
XLk(k\k)  in  the  local  coordinate  frame.  In  the  present  context  (i.e.,  for  negative  time  update) 
*-Lk{k\k)  is  no  longer  the  sufficient  statistic,  and  determination  of  the  optimal  -  i.e.,  the  minimum 
mean  squared  error  (MMSE)  -  estimate  requires  the  knowledge  of  past  data  (i.e,  xff(&  —  1| k  —  1), 
zg(k  —  1)  etc.)  Since  only  the  updated  state  and  covariance  are  stored  our  implementation,  the 
following  sub-optimal  approach  is  used  for  negative  time  updates. 

We  “predict”  the  state  at  tk+i  by  propagating  the  state  equation  backwards  without  process 
noise  i.e., 

x/.Jfc  +  ]  \k)  =  <3>i:XLfc(/:|k)  (25) 

PL,(k  +  l\k)  =  +  (26) 

where  4>fc  is  the  (backwards)  transition  matrix  from  tk  to  tk+i .  Since  the  latest  time  instant  is  tk 
and  not  tk+ 1,  the  desired  updated  state  should  be  at  tk  (in  the  standard  case  it  is  at  tk+ 1).  The 


state  is  updated  directly  using  X-Lk{k  +  l|fc)  as  the  “predicted”  state,  i.e., 

vLk{k+  1)  =  zLk(k  +  1)  -  HxLk(k  +  l|fc)  (27) 

SLk(k+  1)  =  H PLk (k  +  1| k)  H'  +  RLk {k  +  1)  (28) 

WLk(k  +  1)  =  PLk(k  +  l\k)H,SLk(k  +  I)”1  (29) 

x-lk(k\k)  =  xLk(k\k)  +  WLk(k  +  l)vLk(k  +  l)  (30) 

Plk{k\k)  =  pLk(k\k)-wLk{k  +  i)sLk{k  +  i)wLk(k  +  iy  (si) 


The  state  estimate  x.Lk(k\k)  and  covariance  PLk(k\k ))  are  now  updated  to  (k\k)  and  P£,k{k\k) 
respectively.  This  method  has  been  found  to  be  quite  satisfactory  in  all  the  instances  where 
negative  time  updates  were  required. 
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MISSION 

OF 

ROME  LABORATORY 


Mission.  The  mission  of  Rome  Laboratory  is  to  advance  the  science  and 
technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
Rome  Lab: 


a.  Conducts  vigorous  research,  development  and  test  programs  in  all 
applicable  technologies; 

b.  Transitions  technology  to  current  and  future  systems  to  improve 
operational  capability,  readiness,  and  supportability; 

c.  Provides  a  full  range  of  technical  support  to  Air  Force  Materiel 
Command  product  centers  and  other  Air  Force  organizations; 

d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  communications,  command  and  control,  intelligence,  reliability 
science,  electro-magnetic  technology,  photonics,  signal  processing,  and 
computational  science. 


The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  Intelligence,  Signal  Processing, 
Computer  Science  and  Technology,  Electromagnetic  Technology, 
Photonics  and  Reliability  Sciences. 


